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Two Photon Decays of Charmonia from Lattice QCD 

Jozef J. Dudck and Robert G. Edwards 
Jefferson Laboratory MS 12H2, 12000 Jefferson Avenue, Newport News, VA 23606, US^ 

We make the first calculation in lattice QCD of two-photon decays of mesons. Working in the 
charmonium sector, using the LSZ reduction to relate a photon to a sum of hadronic vector eigen- 
states, we compute form-factors in both the space-like and time-like domains for the transitions 
rjc — > 7*7* and XcO 1*1* ■ At the on-shell point we find approximate agreement with experimen- 
tal world- average values. 

PACS numbers; 12.38.Gc, 12.40.Vv, 13.25.Gv, 14.70.Bh 
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The two photon widths of charmonia have attracted 
the attention of a range of a theoretical and experi- 
mental techniques. Within perturbative QCD the two- 
photon branching fraction is believed to give access to 
the strong coupling constant at the charmonium scale 
though cancellation of non-perturbative factors, while in 
quark models and the effective field theory NRQCD, two- 
photon widths have been proposed as a sensitive test 
of the corrections to the non-relativistic approximation. 
Experimental measurements are diverse, coming from on 
the one hand two-photon fusion at e^e~ machines with 
subsequent reconstruction of the charmonium in light 
hadrons, and, on the other, pp annihilation to charmo- 
nium with decay to real 77 pairs being detected. Im- 
provements in extracted branching fractions from both 
methods will soon arrive via more accurate measurements 
from Belle, BaBar and CLEO-c of exclusive hadronic rjc 
decays. 

Despite the considerable theoretical and experimental 
interest, there has not yet been an ab-initio estimation 
of two photon widths in charmonia. In this paper we ad- 
dress this issue using a novel application of lattice QCD. 

Conventionally lattice QCD calculations involve the 
evaluation of the matrix element of an operator between 
hadron states. The hadrons are induced using interpolat- 
ing fields (some combination of quark and gluon fields) 
with the appropriate hadron quantum numbers. However 
these fields will typically have overlap with many excited 
states which have the same external quantum numbers 
so that in order to isolate the ground states, the interpo- 
lating fields are taken out to Euclidean times far from the 



operator (and each other). Such an approach will clearly 
not work if the initial or final state contains no hadron 
as is the case for two-photon decays since the photon is 
not an eigenstate of QCD. In this case we can adopt a 
slightly more sophisticated approach, using the formal 
relationship between the S-matrix and field-theoretic TV- 
point functions and the accurate perturbative expansion 
in the photon-quark coupling to express the photon as a 
superposition of QCD eigenstates. 

The concept and realization is lucidly presented by Ji 
and Jung in where they consider the hadronic struc- 
ture of the photon using lattice three-point functions. We 
will outline the method as applied to meson two-photon 
decays. 

The amplitude for the two-photon decay of a meson M 
can be expressed in terms of a photon two-point function 
in Minkoswki space by means of the LSZ reduction 

(7('?i,Ai)7(g2,A2)|M(p)) = - lim e;(gi, Ai)e:((72, A2) 
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d^xd^ye"i'^-y+"i'^-''{0\T{A^'{y)A%x)}\M{p)), 



up to photon renormalisation factors. The explicit pho- 
ton fields prevent direct computation of this quantity in 
lattice QCD, however we can utilize the perturbative ex- 
pansion of the photon-quark coupling to approximately 
integrate them out (the path-integral over gluon fields is 
suppressed): / VAVi^V^^e'-^'^'^^'^'^^'i'^'l-'^ At'{y)A'' [x) = 



j PAPV;DV'e'^°[^''^''^](... + 4 j dt^zdt^w [4,{z)Yi}j{z)Ap{z)\ [4,{w)i'' iP{w)A^{w)\ + . . .) A^' {y) A" {x) , 

The integration over the photon field can be carried out by Wick contracting the fields into propagator products, so 
that, neglecting disconected pieces, {i{qi,\i)i{q2T\2)\M{p)) = 

(-e^) hm e*^{qu\i)el{q2A2)q?q2 [ d^x d^w d'' z e^-'^ D'^P{0, z)D'''^ {x, w){0\T{jpiz)Mw)}\Mip)) . 
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The photon propagator can be written D'"^(0, z) 
the integral. 



(27r)'i k'^+it 1 cancelling the inverse propagators outside 



As explained in , the resulting expression can be ana- 
lytically continued from Minkowski to Euclidean space- 
time provided the photon virtualities, Q\ = — 
Q2 = |g2p — are not sufficiently timelike that 



they can produce on-shell hadrons. In charmoniumjl 1| 



r 2 eA'('?i,'^i)ei^(92, A2) 
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this limits us to > 



Using suitable a 



QCD interpolating field to produce M and reversing 
the operator time-ordering for convenience we have 
(A'/(p)|7((7i,Ai)7(g2, A2)) = 



d^xe-^''-%M{x,tf) [ d^ye^^-^yf{y,t)f{Q,U)]\Q) (1) 
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It is clear from the previous discussion that obtaining 
two-photon widths is a natural extension to the study 
of radiative transitions carried out in Q - there we 
computed three-point functions involving vector currents 
with the source {ti) and sink {tf) positions fixed and var- 
ied the vector current "insertion" [t] across the temporal 
direction to plot out a plateau. For two-photon widths 
we repeat this but with a varying sink (or in the case of 
eqn. source) position which will be integrated over 
with an exponential weighting. 

Our calculation was performed in the quenched trun- 
cation of QCD using the Chroma software system. We 
employed a 24"^ x 48 lattice generated using a Wilson 
gauge action with coupling (3 = 6.5 and lattice spacing 
a « 0.047 fm (determined from the static quark potential 
in [3). Quark propagators were computed using a non- 
perturbatively improved Clover action with Dirich- 
let boundaries in the temporal direction, with the charm 
quark mass set using the spin-averaged 77c, J/V" mass, ac- 
curate to within 3%. Since we adopted only smeared 
local fcrmion bilinears as meson interpolating fields we 
had access to three-point functions involving 77c and XcO 
mesons. In typical radii of charmonium states were 
extracted and foimd to be much smaller than the ~ 1.2 
fm box used here. We adopt the conserved ( "point-split" ) 
vector current at the insertion - our experience in radia- 
tive transitions is that including the improvement term^] 
has an effect below 5% near = 0. 

We applied two different methods to calculate . The 
first (using 174 configurations) was to place the meson 
state at a fixed sink position = 37. As in the sink 
was used as a sequential source for a backward propa- 
gator inversion, meaning that its properties were fixed 
for each computation while we were able, without fur- 
ther cost, to vary the direction and momentum of the 
insertion and the direction of the source field. We then 
computed with all possible source positions, ti, which, 
while costly in computing time, allowed us to freely vary 



the value of wi and hence Q\ and in addition to view the 
subsequent integrand. In figure ^a) we display the inte- 
grand for "insertion" positions t ^ 4, 16, 32, pf ~ (000) 
and qi = (100) with an rjc at the sink. 

It is clear that provided the insertion is not placed 
too close to the Dirichlet wall (i.e. i ^ 7) we will 
be able to capture the full integral by summing time 
slices, ti. The integral as a function of insertion posi- 
tion t is shown in figure ^fb) for a selection of Q2 with 
Qf = {uji = \qi\) where we observe plateaus with 
the deviation from plateau behavior at larger t com- 
ing from excited 77^ contributions, both of which are fit- 
ted simultaneously. Extracting the plateau values for 
a range of Qi (which we are free to choose continu- 
ously) and Q2 (which is fixed for a given set of wi, g2,p), 
we find the dependence displayed in figure El We 
plot dimensionless F defined by {rid'yiqi, Xi)j{q2, ^2)) = 
2{lefm-}F{QlQl)e^^p^e'^e^q^q^, where the on-shell 
decay width is r(77c — > 77) = TTaljy^^mnJF{0,0)\'^. 

In order to more clearly express the voluminous data 
in figure El we adopt a simple one-pole parameterization 
to fit the data for each value of Ql, 



F{QlQl) 



FjQlO) 
1 ^ Ql 



(2) 



The curves in figure El are fits of this form. In the insets 
of figurcElwe display the fit parameters F{Qi, 0), /i(Qj). 
F{Q 1,0) is itself then fitted with a one-pole form with 
F(0, 0) = 0.168(4), /t = 3.040(20) GeV, which is compati- 
ble with ij,{0) = 3.008(16) GeV indicating the reasonable- 
ness of the one-pole parameterization in this region. 
The both-photons on-shell point is three standard devi- 
ations smaller than the value extracted from the PDG0 
decay width = 0.274(43). 

The error quoted on our result is statistical only and 
must be augmented by an uncertainty due to scaling from 
our fixed lattice spacing to the continuum and one related 
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FIG. 1; (a) Integrand in at three values of vector cur- 

rent insertion time {t = 4, 16, 32) with pseudoscalar sequential 
source at sink position t/ = 37. (b) Pseudoscalar two- photon 
form- factor as a function of time slice, t, from First six 

time slices ghosted out due to the Dirichlet wall truncating 
the integral. Constant plus single exponential fits shown in 
orange. 



to the lack of light-quark loops within the quenched trun- 
cation. Since we compute at only one lattice spacing we 
cannot accurately determine a scaling uncertainty, how- 
ever we expect the non-perturbatively tuned Clover ac- 
tion to have small 0{a) scaling errors. On the other hand 
it is possible that we have non-negligible 0(mca) scaling 
which we conservatively estimate to be at the 15% level. 

In the NRQCD factorized approach the decay ampli- 
tude is proportional to the rjc wavefunction at the ori- 
gin with relativistic corrections. In the quenched the- 
ory utilized here, the incorrect running of the coupling 
down to short distances with the scale set at long dis- 
tances, does typically lead to a depletion of the wavefunc- 
tion at the origin. On the lattice used here we find[l^ 
/,/, = 386(6) MeV in comparison to the experimental 
411(7) MeV, however we cannot determine to what de- 
gree the quenched depletion is offset by 0{mca) scaling 
effects, and as such we assign an estimated 10% error due 
to quenching. 

The similarity of /i(0) to the J/V' mass in our calcula- 
tion (3006(5) MeV) suggests an approximate vector me- 




FIG. 2: two-photon form-factor, F{Q\,Q2)- Points are 
lattice QCD data, fits are monopole forms as described in the 
text. Lower left shows fitted amplitudes at Qi = and lower 
right the fitted pole masses /i(Qf). 



son dominance description of the form-factor. This is 
at first sight somewhat surprising given that VMD was 
demonstrated not to hold for the r]c charge "form-factor" 
in 0. This was explained as being due to the presence 
of very many closely spaced vector poles (the J/'i/', i/j' ■ ■ ■) 
which must be summed - in contrast to the light-quark 
sector where the p(1450) is rather distant in compari- 
son to the p{770) and where VMD tends to be a good 
approximation. 

What we suspect is happening here is that while the ex- 
cited ip poles are not negligible by their remoteness, they 
are almost negligible by their small residues. Consider 
the case = and the time ordering in which photon 
1 is emitted from the r]c first, then the time-dependent 
perturbation theory amplitude is proportional to 



E 

N 



(^?cb-;.(o)|V>W)(V>W|>(o)|o) 



(3) 



The numerator is essentially the residue of the A^th ex- 
cited '0 pole and is seen to be proportional to the prod- 
uct of the Ml transition matrix element between rjc and 
-0'^^^ and the decay constant of the ip'-^\ The decay 
constants fall slowly with increasing N as observed from 
the experimental e~^e~ widths, while the Ml transition 
amplitudes are expected to fall rather rapidly - they are 
"hindered" transitions that are proportional to the over- 
lap of orthogonal wavefunctions with a small correction 



for recoil. Because of this only the J/^/i term in the sum 
has a considerable residue (since the wavefunction over- 
lap with the is close to one) and we observe something 
like VMD. A hint of this behavior was observed in the 
L3 experiment , and a model realization demonstrated 
in 0. 

The fact that ^ varies with Ql suggests that as this 
photon goes off-shell the higher ip resonances get a larger 
residue (through a larger recoil correction to the hindered 
transition) and contribute more to the sum over poles 
moving the effective pole position, fJ.{Ql)- 

The method of computing three-point functions with 
all possible source positions {U = 0...tf) is extremely 
costly in computing resources. With the penalty of los- 
ing the ability to freely vary Qi we can reduce the com- 
puting time by a factor of 0{Lt) by putting the meson 
interpolating field at the source and using 



LOlt 



(4) 



in the sequential source for the backward propagator in- 
version. It is then necessary to fix wi and qi in advance 
and one cannot view the integrand since the integra- 
tion is being performed "on-the-fly" within the sequential 
source. We computed 300 configurations with = in 
this way and the results for 77c and Xco are displayed in 
figure 13 Note that the one-pole fit to the rjc data yields 
F(0,0) = 0.160(8) and fi = 2.899(71), values that are 
in agreement with those extracted from the method de- 
scribed above. 

We define the XcO two-photon form-factor by 



(Xco|7(9i,Ai)7(g2,A2)) -2(|e)2 

I £1 • e2(?l • 92 



xm-^^G{QlQl){e 



£2 • giei • qi 



(5) 



The PDG|3 average of experimental measurements 
for the two-photon width corresponds to |G(0,0)| = 
0.150(13) with which we are in excellent agreement - our 
single-pole fit to lattice data returns |G(0, 0)| = 0.146(16) 
and At = 2.976(81) GeV. 

The logic used to justify the approximate J/i/j VMD 
observed for the ijc does not so obviously apply to the 
XcO- The radiative transitions to virtual vector meson 
states are now El and are not "hindered" for the excited 
states to the same extent as Ml (This can be seen in, for 



(Thi 

■m) 



example, the quark model study [lOj ) . Since the residues 
of excited state poles are not falling so quickly they can 
contribute to the sum and wc will not necessarily see 
dominance of the J/tp. 

Unfortunately in this case the two nearest poles, the 
-0(3686) and the 0(3770) can conspire to cancel each 
other's effect. Using a combination of experimental data 
and quark model predictions one estimates that the 
residues of these two poles are approximately equal (they 
have value of about 1/3 the J/tp residue). If these 
residues happen to have opposite sign then the J/ip pole 




FIG. 3: (a) 77c — > 77* amplitude, (b) XcO 77* amplitude. 
Fits are monopole forms as described in the text. 



will dominate with small contributions from the more 
highly excited ip resonances. Our data is insufhciently 
precise to distinguish any deviation from a one-pole be- 
havior. 

In summary we have demonstated the feasibility of 
using a suitable sum over lattice timeslices to sim- 
ulate a photon in an external state, using the phe- 
nomenologically interesting case of the two-photon de- 
cays of charmonia. Our results r(r7c 77) = 
2.65(26)stat(80).cai.(53)quen.keV, r(xco ^ 77) = 

sys 

2.41(58)stat(72)scai.(48)quon. kcV, are in reasonable agree- 

sys 

mcnt with experiment. We believe them to be systemat- 
ics dominated, in particular owing to the neglect of light 
quark loops and discretisation effects related to the heavy 
quark mass. The systematic error can be reduced with 
further computations, which arc now warranted given the 
feasibility demonstrated herein. This lattice technique 
has the advantage over non-rclativistic models that it can 
be applied to the light-quark sector without fundamental 
change. 
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publish or reproduce this manuscript for U.S. Govern- 
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